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' Abstract 

(N. In this literature, we carefully investigate the structure of single- and multi-frequency imag- 

ing functions, that are usually employed in inverse scattering problems. Based on patterns of 
the singular vectors of the Multi-Static Response (MSR) matrix, we establish a relationship 
• between imaging functions and the Bessel function. This relationship indicates certain prop- 

erties of imaging functions and the reason behind enhancement in the imaging performance 
by multiple frequencies. Several numerical simulations with a large amount of noisy data are 
performed in order to support our investigation. 

i ■ 1 Introduction 
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One of the main objective of the inverse scattering problem is to identify the characteristics of 
unknown targets from measured scattered field or far-field pattern. In research fields such as 
physics, medical science, and materials engineering, this is an interesting and important problem. 
Related works can be found in [3J [El [HI H2J 1131 HH] and references therein. In order to solve 
this problem, various algorithms for finding the locations and/or shapes of targets have been 
accordingly developed. 

In many research studies [SJ [8j [TTJ [14j [18j [27] , the shape reconstruction method is based on 
Newton-type iterative algorithms. However, for a successful shape reconstruction using these 
algorithms, the iterative procedure must begin with a good initial guess that is close to the unknown 
target because it highly depends on the initial guess; for more details, refer to [T71 |2"T] . 

To finding a good initial guess, alternative non-iterative reconstruction algorithms have been 
developed, such at the Multiple Signal Classification (MUSIC)-type algorithm [6j [24j [26] , linear 
sampling method [§J HH] , topological derivative strategy [H [7J [THl HH H3J , and the multi- frequency 
based algorithm [31 [SJ [101 dU [25] . Among them, although the multi- frequency based algorithm 
has exhibited potential as a non-iterative imaging technique, a mathematical identification of its 
5^ 1 structure needs to be performed for its heuristical applications, which is the motivation behind. 

In this paper, by intensively analyzing the structure of single- and multi- frequency imaging 
functions, we discover some properties of imaging functions and confirm the reason behind the 
enhancement in the imaging performance by applying multiple frequencies. In recent work [5J, this 
fact was verified by the Statistical Hypothesis Testing but our approach is to find a relationship 
between imaging functions and Bessel functions of the first kind of the integer order. 

This paper is organized as follows. In section [H we briefly review the two-dimensional direct 
scattering problem, and an asymptotic expansion formula for far-field patterns, and introduce 
the imaging function introduced in |20j . In section [3J we analyze the single- and multi- frequency 
based imaging functions and discuss their properties. In section 21 we present several numerical 
experiments and discuss the effectiveness, robustness, and limitation of imaging functions. Finally, 
a brief conclusion is given in section [5] 
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2 Review on imaging function 



In this section, we survey the two-dimensional direct scattering problem and an imaging algorithm. 
A more detailed dcscussion can be found in [201 EH 125] . 

2.1 Direct scattering problem and asymptotic expansion formula 

Let S m be a homogeneous inclusion with a small diameter p in the two-dimensional space R 2 . 
Throughout this paper, we assume that every E m is expressed as 

where r m and p denote the location and size of E m , respectively. Here, T> m is a simple connected 
smooth domain containing the origin. 

Let £o and po respectively denote the dielectric permittivity and magnetic permeability of 
IR 2 . Similarly, we let e m and p m be those of E m . For simplicity, let E be the collection of E m , 
m = 1, 2, ■ • • , M, and we define the following pieccwisc constants: 

, , _ J £ m for r G E m . . _ J p m for re E m 

£[r >-\ e for r£R 2 \E and W> ~ \ Mo for r £ R 2 \E. 

At a given frequency w, let u tot (r, u>;l) be the time-harmonic total field that satisfies the 
Hclmholtz equation 

V- (-^Vu M (r,u-,l)j +uj 2 e(r)u tot (r,uj;l) =0 (1) 

with transmission conditions 

Let u ino (r, lo;1) be the solution of (JlJ without E. In this paper, we consider the following 
plane-wave illumination: for a vector d/ £ C 1 , u inc (r,u;;Z) = exp(jwd/ ■ r). Here, C 1 denotes a 
two-dimensional unit circle. 

Generally, the total field u tot can be divided into the incident field u inc and the unknown 
scattered field u acat , which satisfies the Sommcrfcld radiation condition 

lim y/\r\ ( — ff^li^iil _ jfc u scat (r,w; Z) ) = 0, k = u}y/e fM> 
|r|^oo \ o|r| j 

uniformly in all directions f = A . As given in [5] , u scat can be written as the following asymptotic 
expansion formula in terms of p 

M 



w 3cat (r,w; I) = p 2 ^ ( Vu inc (r m ,cj;/) • T(r m ) • V$(r m ,r, u) 



m—1 

+ u 2 (e- e )u mc (r m , u; l)<f>(r m , r, w)J + o(p 2 ), (2) 

where o(p 2 ) is uniform in r m £ E m and d; £ <£ . Here $(r r „,r, w) is the two-dimensional time 
harmonic Green function (or fundamental solution to Hclmholtz equation) 

$(r m ,r,w) = -^ ^-ffo( fc o|r m - r|), 

where Hq is the Hankcl function of order zero and of the first kind. 
The far-field pattern is defined as function F(r, d) that satisfies 

exp(jfc |r|) / 1 \ 

u Bcat (r,u;;Z) = = F(r,d) + o — = (3) 



oo uniformly on f = A. 



2 



2.2 Introduction to imaging function 

The imaging algorithm introduced in |20j used the structure of a singular vector of the Multi-Static 
Response (MSR) matrix M = (F pq ) = (F(r p ,d q ))^ q=1 . Note that by combining ©, and the 
asymptotic behavior of the Hankcl function, the far-field pattern F pq with observation number p 
and incident wave number q can be represented as the asymptotic expansion formula (see [5] for 
instance) 



M 



F(r p ,d q )^p 2 



V V £ oMo 



4^ 



- r p ■ T(r m ) • d g exp jfco(d 9 - r p ) • r 



(4) 



Here T(r rn ) is a 2 x 2 symmetric matrix: 

T(r m ) : 



2^o 



Mm + MO 



area(£> m )I 2 , 



where I„ denotes the n x n identity matrix. For the sake of simplicity, we eliminate the constant 
4y/k J * n ®' Then, the incident and observation direction configurations are kept same, i.e., for 
each r p = — d p , the pq-th element of the MSR matrix M is given by 

M r 



=— d„ 



p 2 E 



2mo 



\/£0M0 Mm + MO 



area(£> m )d p • dq 



exp ifco(d p + d g ) • r Tl 



where area(2? m ) denotes the area of T> m . 

Based on the above representation of F pq , we introduce a vector D(r;u;) e C Wxl as 



/ ei exp(jfc di • r) ^ 
e 2 exp(jfc d 2 • r) 
D(r;w):= . 

\ e w cxp(jfc dAr ■ r) / 

Then M can be decomposed as follows: 



where e p = (1, d p 



(5) 



M 



so 



M= V D(r m ;w) ^JIb~ D(r m ;w) T . 

m=l V T ( r ™) / 

This decomposition leads us to introduce an imaging algorithm as follows. First, let us perform 
the Singular Value Decomposition (SVD) as follows: 



M 



(6) 



where U m and V m are the left and right singular vectors, respectively. Then, by defining a 
normalized unit vector 

D(r;w) 



W(r;W) -|D(r; W )r 
we can introduce an imaging function as follows (see [51 l2"01 |2"T1 12"5] ) 

M 



J2 (w(r;«) • U m (w)) (w(r;w) • V m (u>)) 



(7) 



(8) 
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3 Structure analysis of imaging function 

3.1 Structure of imaging function 



We now determine! the structure of imaging function ([5)1. For this purpose, we recall some useful 
statements. 

Lemma 3.1 ([5]). A relation A ~ B means that there exists a constant C such that A = CB. 
Then, for vectors U m and V m in (0) and W(r;w) in the following relationship holds 

W(r m ;u;) ~ U m (uj) and W(r m ;oj) ~ V m (u)). 

Lemma 3.2 ([15]). Let d,r G R 2 , ant! w > 0; then 



/ exp(jwd ■ r)dS(d) = 27rJ (w|r|), 



where J u (x) denotes the Bessel function of order v of the first kind. 

Subsequently, we can explore the structure of ([5]) as follows 

Theorem 3.3. If the number of incident and observation directions is sufficiently large, then the 
imaging function can be represented as follows: 



M 



(9) 



Proof. By hypothesis, we assume that N is sufficiently large. For simplicity, we set Ad p := 
| dp — d p _i| for p = 2, 3, ■■• ,N, and Adi := |di — djv|. Ad, is defined analogously. Then, 
applying Lemmas 13.11 and 13.21 yields 



r; w = 



M 



m—1 

M 



J2 W(r; u) ■ UM W(r; w) • V m (w) 



N 



Y Y cx p0' wd p • ( r m - r )) 



m—1 \p— 1 

M 



Adp 
2tt 



AT 



^exp(jwd g ■ (r m - r)) 



Ad, 

2tt 



(10) 



1 

4^ 



V] / exp(jo;d • (r m - r))dS(d) 



E J o( 



w rv 



This completes the proof. 



□ 



Note that Jq(x) has the maximum value 1 at x — 0. This is the reason why the map of W(r; w) 
plots magnitude 1 at r = r m € £ m . Moreover, due to the oscillating property of Jo(x), Theorem 
3.31 indicates why imaging function (jSJ plots unexpected replicas, as shown in Figure [TJ 



3.2 Reason behind enhancement in the imaging performance by apply- 
ing multiple frequencies 



According to the Theorem 13.31 the oscillating pattern of the Bessel function must be reduced or 
eliminated in order to improve the imaging performance. One way to do so is to apply the high- 
frequency lj = +oo in theory. Another way is to apply several frequencies to the imaging function 
dHI) as follows: 



(r;S) := 



1 



]Tw(r; Ws ) 



s=l 



M 



s—1 m—1 



Y W ( r ' ^ ■ U '»( w s) W(r; u) s ) 



(11) 
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Several researches in [31 [SJ HOI HI] have confirmed on the basis of Statistical Hypothesis Testing 
and numerical experiments that the multi-frequency imaging function (jllj) is an improved version 
of the single- frequency version (O. The reason for this is discussed as follows. 

Theorem 3.4. If u>g and the number of incident and observation directions are sufficiently large 
enough, then the structure of the imaging function 111]) is 



W(r; S) 



J 2 (w s |r m - r|) + J\ (uj s \r m - r|) 



uj s - u>i 



Proof. According to ((TU]), we can observe that 

s 



(r;5) 



]T J$(us\r m -r|) 



s=l 



uj s - wi 



Jo(wi|r m - r|) + Jf(uii\r m - r|) 



J${u)\r m - r\)du 



(12) 



Using this, wc apply an indefinite integral formula of the Bcssel function (see HE1 page 35]): 

J 2 (x) = x( J$(x) + Jf(x)) + f Jl{x)dx 



in addition to a change of variable w|r m — r| = x. This yields 



Jo(w|r m - r|)cfcj = | 1 / Jl(x)dx 



r m - r 



Wl|r m — r| 



=cj s I Jo(wg|r m - r|) + Ji(u} S \r m - r|) ) - cjW Jo(a;i|r m - r|) + J 2 (wi|r m - r|) 



/ Ji(w|r m - r|)dc 



Now, we consider the upper bound of 



A(|r m - r|,wi,ws) 



Jf(u}\r m - r\)dw. 



Note that since Ji(0) = 0, let us assume that |r m — r| 7^ and < uis\r m — r| <C v2. Then 
applying asymptotic behavior 

rn 1 ( ^ " 

J v [x) 



r(i/ + i) V2 

and boundedness property ,7„(a;|r m — r|) < yields 



us ^ Z^s 1 rus\r m -r\ J 1 ( x ) 
jf (w|r m — r|)dw < — p / Ji(w|r m — r|)dw w — p / r- rcfe 



a/2 J Wl 
1 



2\/2|r TO - r| ,/ Wl |r ra -r| 



"s|r m -r| 



\/2 y wi |p m -r| l r m 

xd^ ( ^ )2 "> )2 |r m -rl 



4\/2 
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Figure 1: (a) 2-D plot of © and (b) 2-D plot of (JTTJ when m = 1 and r m = 0. 



Now, assume that cos satisfies 



V2 



Then since 



we can obtain 



ws|r m -r| > V2 i.e., |r m - r| > — > 0. 



J\{x)dx = —Jq(x), 



Ji(u\r m - r\)du> 



w s |r m -r| 



Wi|r m — r| 
f 



dx 



ui\r m — r| 



< 



|r m - r| 
2 



J (wi|r m - r|) - J {io s \r m - r| 
< V2w s . 



Therefore, the term A(|r m — r|,o;i,Wsr) can be disregarded because 



uj s - UJl 



J (w s |r m - r|) + J l {u s \r m - r|) 



J 2 (wi|r m - r|) + Jl(ui\r m - r|) 



uj s - u)i 

and A(|r m — r\,u>i,u)s) "C O(ws). Hence we can obtain ([12")) . This completes the proof. 



0(w s ) 



□ 



The graph for l|12p is shown in Figure [T] This graph shows that (fTT|> yields better images 
owing to less oscillation than Q docs. This result indicates why a multi-frequency based imaging 
function offers better results than a single-frequency based one. 
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Figure 2: Maps of W(r; 10). 



4 Numerical experiments and discussions 

In this section, we describe the numerical experiments we conducted to validate our analysis. For 
this purpose, we choose a set of three different small disks E m . The common radii p of S m are set 
to 0.1, and parameters Eq and [1q are chosen as 1. Locations r m of S m are selected as 

ri = (0.4, 0), r 2 = (-0.6, 0.3), and r 3 = (0.1, -0.5). 

For a given wavelength A s , each frequency is selected as u s = for s = 1, 2, • • • ,5. Note that all 
the wavclcngthcs A s arc uniformly distributed in the interval [Ai,As]. The observation directions 
dp G S 1 are selected as 

/ 2irp , 2irp\ 

d„ = cos , sm tor p = 1, 2, • • • , N. 

\ N N J y ' ' 

The incident directions d q G S 1 are selected analogously. 

In all the examples, a white Gaussian noise with 10 dB signal-to-noise ratio (SNR) is added 
to the unperturbed data in order to exhibit the robustness of the proposed algorithm via the 
MATLAB command awgn. In order to obtain the number of nonzero singular values M for each 
frequency oj s , a 0.01-threshold scheme is adopted (see [231 US]). The search domain f2 is selected 
as a square D. = [—1, 1] x [—1, 1]. 

Figure^] shows the map of W(r; 10) via the MSR matrix M for N = 20 and S = 10 and different 
frequencies with Ai = 0.5 and As = 0.3. On the left-hand side of Figure[2l we set the same material 
properties e rn = 5 and /i m = 5, m — 1, 2, 3. As expected, locations of E m can be clearly identified. 
On the right-hand side of Figure we set different material properties e\ = fix = 5, £2 = A*2 = 2, 
and £3 = /13 = 7. Note that due to the small values of £2 and 112, the map of W(r; 10) plots a small 
magnitude at 2 but the locations of all S m are well identified. 

Figure [3] shows the influence of the number of applied frequencies S. As we discussed in section 
[3l increasing S yields a more accurate image. Note that applying an infinite number of S would 
yield good results in theory, but in this experiment, S = 10 is sufficient for obtaining a good result. 

Based on the recent work |20j . the proposed algorithm can be applied to imaging of extended, 
crack-like electromagnetic inclusion(s) T with a supporting curve 7 and a small thickness h. How- 
ever, note that even a sufficiently large number of N and S applied to obtain a good image of 
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Figure 3: Maps of W(r; S) with (a) S 



= 1 (b) S 



= 3 (c) S = 7 and (d) 5 



20. 




Figure 4: Maps of W(r; S). 



a complex-shaped thin inclusion using the proposed algorithm can occasionally yield poor results 
(see Figure HJ. Furthermore, note that the elements of M are expressed as written by 



F(d p ,d q )~Yl 



1 M, , 

d„ -t 



£^ + 2 ( - - — ) d p ■ t(r m )d, • t(r m ) 

+2[ — - \ )d p • n(r m )d q • n(r m ) 



exp jfc (d p + dq) • r ri 



Therefore, e p in ([5]) must be a linear combination of a unit tangential vector t(r m ) and a normal 
vector n(r m ) at r G 7. If we have a priori information of t(r m ) and n(r m ), we can obtain a good 
result. However, because this is not the case, it is difficult to obtain a good result. This is further 
explained in detail in Section 4.3.1]. 



5 Conclusion 

Using an integral representation formula and an indefinite integral of the Bessel function, we 
determined the structure of single and multiple electromagnetic imaging functions. Because of the 
oscillation aspect of the Bessel function, we confirmed the reason behind the improved imaging 
performance by successfully applying high and multiple frequencies. 

Recent works [20j [2lJ [25] have shown that the proposed algorithm can be applied to limited- 
view inverse scattering problems. In order to determine the structure of imaging functions, the 
integration in Lemma |3~21 on the subset of a unit circle must be evaluated; however, this evaluation 
is very difficult to perform. Therefore, identifying the imaging function stucture in the limited- 
view problem will prove to be an interesting research topic. Moreover, in the imaging of crack- 
like inclusions, estimating unit tangential and normal vectors on such an inclusion and yielding 
relatively good results will be an interesting work. 

Finally, we considered the imaging function for penetrable electromagnetic inclusions but it will 
be applied to the perfectly conducting inclusion(s) directly. Extension to the perfectly conducting 
target will be a forthcoming work. We further believe that the proposed strategy can be extended 
to a three-dimensional problem. 
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